#!/usr/bin/env python3
"""
Geometric Pi Calculator

This script implements various geometric and integration-based methods for calculating π.
These methods visualize the mathematical concepts behind π and demonstrate the
relationship between π and various geometric shapes and integrals.

Methods implemented:
1. Circle Approximation - using random point sampling with visualization
2. Integration Methods - using numerical integration techniques
3. Polygon Perimeter Method - with incremental visualization
4. Historical Approximations - implementations of historical π approximations

The script also includes visualization capabilities to help understand
the geometric intuition behind π.
"""

import math
import random
import time
import numpy as np
import sys
from matplotlib import pyplot as plt
from matplotlib.patches import Circle, Polygon
from decimal import Decimal, getcontext


def circle_approximation_method(points=1000, visualize=False):
    """
    Approximate π using a circle inscribed in a square with point sampling.
    
    This method visualizes the Monte Carlo approach by actually plotting
    the points in a coordinate system.
    
    Args:
        points: Number of random points to generate
        visualize: Whether to create a visualization
        
    Returns:
        Approximation of π and optionally displays a plot
    """
    inside_circle = 0
    x_inside, y_inside = [], []
    x_outside, y_outside = [], []
    
    # Generate random points and check if they're inside the unit circle
    for _ in range(points):
        x = random.uniform(-1, 1)
        y = random.uniform(-1, 1)
        
        if x**2 + y**2 <= 1:
            inside_circle += 1
            x_inside.append(x)
            y_inside.append(y)
        else:
            x_outside.append(x)
            y_outside.append(y)
    
    # Calculate π: area of circle = π*r², area of square = 4*r²
    # So π = 4 * (points inside circle / total points)
    pi_approx = 4 * inside_circle / points
    
    # Create visualization if requested
    if visualize:
        plt.figure(figsize=(10, 10))
        
        # Plot the unit circle
        circle = plt.Circle((0, 0), 1, fill=False, color='black', linewidth=2)
        plt.gca().add_patch(circle)
        
        # Plot the points
        plt.scatter(x_inside, y_inside, color='blue', alpha=0.5, label='Inside Circle')
        plt.scatter(x_outside, y_outside, color='red', alpha=0.5, label='Outside Circle')
        
        # Draw the square
        plt.axhline(y=-1, color='black', linestyle='-', linewidth=2)
        plt.axhline(y=1, color='black', linestyle='-', linewidth=2)
        plt.axvline(x=-1, color='black', linestyle='-', linewidth=2)
        plt.axvline(x=1, color='black', linestyle='-', linewidth=2)
        
        # Add labels and title
        plt.title(f'Monte Carlo Estimation of π: {pi_approx:.6f} (Error: {abs(pi_approx - math.pi):.6f})')
        plt.xlabel('x')
        plt.ylabel('y')
        plt.legend()
        plt.axis('equal')
        plt.grid(True)
        plt.savefig('monte_carlo_pi.png')
        plt.close()
        print(f"Created visualization 'monte_carlo_pi.png'")
    
    return pi_approx


def integration_methods(method='midpoint', intervals=1000, visualize=False):
    """
    Calculate π using numerical integration methods.
    
    Uses the fact that ∫₀¹ 4/(1+x²) dx = π
    
    Args:
        method: Integration method ('midpoint', 'trapezoid', or 'simpson')
        intervals: Number of intervals for numerical integration
        visualize: Whether to create a visualization
        
    Returns:
        Approximation of π
    """
    def f(x):
        return 4 / (1 + x**2)
    
    a, b = 0, 1  # Integration range
    h = (b - a) / intervals
    
    if method == 'midpoint':
        # Midpoint method
        result = 0
        for i in range(intervals):
            x = a + h * (i + 0.5)
            result += f(x)
        pi_approx = h * result
    
    elif method == 'trapezoid':
        # Trapezoid method
        result = f(a)/2 + f(b)/2
        for i in range(1, intervals):
            x = a + h * i
            result += f(x)
        pi_approx = h * result
    
    elif method == 'simpson':
        # Simpson's rule
        result = f(a) + f(b)
        for i in range(1, intervals, 2):
            x = a + h * i
            result += 4 * f(x)
        for i in range(2, intervals-1, 2):
            x = a + h * i
            result += 2 * f(x)
        pi_approx = h * result / 3
    
    else:
        raise ValueError(f"Unknown integration method: {method}")
    
    # Create visualization if requested
    if visualize:
        x = np.linspace(0, 1, 1000)
        y = [f(xi) for xi in x]
        
        plt.figure(figsize=(10, 6))
        plt.plot(x, y, 'b-', linewidth=2, label='f(x) = 4/(1+x²)')
        
        # Shade the area under the curve
        plt.fill_between(x, y, color='skyblue', alpha=0.4)
        
        # Add labels and title
        plt.title(f'{method.capitalize()} Integration: π ≈ {pi_approx:.6f} (Error: {abs(pi_approx - math.pi):.6f})')
        plt.xlabel('x')
        plt.ylabel('f(x)')
        plt.grid(True)
        plt.legend()
        plt.savefig(f'{method}_integration_pi.png')
        plt.close()
        print(f"Created visualization '{method}_integration_pi.png'")
    
    return pi_approx


def polygon_perimeter_method(n_sides=6, iterations=10, visualize=False):
    """
    Calculate π using polygon perimeter approximation with visualization.
    
    This method starts with a regular hexagon inscribed in a unit circle
    and progressively doubles the number of sides to approach π.
    
    Args:
        n_sides: Initial number of sides (default is 6)
        iterations: Number of iterations to double the sides
        visualize: Whether to create visualization
        
    Returns:
        Final approximation of π
    """
    # Start with a regular polygon
    current_sides = n_sides
    
    # For a circle with radius 1, the perimeter of the inscribed polygon
    # with n sides can be calculated using the sine function
    side_length = 2 * math.sin(math.pi / current_sides)
    perimeter = current_sides * side_length
    
    # Array to store results for visualization
    results = [(current_sides, perimeter)]
    
    # Iteratively double the number of sides
    for _ in range(iterations):
        current_sides *= 2
        
        # New side length when doubling the sides
        side_length = 2 * math.sin(math.pi / current_sides)
        perimeter = current_sides * side_length
        results.append((current_sides, perimeter))
    
    # Create visualization if requested
    if visualize:
        plt.figure(figsize=(15, 10))
        
        # Plot perimeter vs. sides
        sides = [r[0] for r in results]
        perimeters = [r[1] for r in results]
        
        plt.subplot(2, 2, 1)
        plt.plot(sides, perimeters, 'bo-')
        plt.axhline(y=2*math.pi, color='r', linestyle='--', label='π = {:.6f}'.format(math.pi))
        plt.xlabel('Number of Sides')
        plt.ylabel('Perimeter')
        plt.title('Perimeter vs. Number of Sides')
        plt.xscale('log')
        plt.grid(True)
        plt.legend()
        
        # Plot error vs. sides
        errors = [abs(p - 2*math.pi) for p in perimeters]
        
        plt.subplot(2, 2, 2)
        plt.plot(sides, errors, 'go-')
        plt.xlabel('Number of Sides')
        plt.ylabel('Error')
        plt.title('Error vs. Number of Sides')
        plt.xscale('log')
        plt.yscale('log')
        plt.grid(True)
        
        # Plot visualization of polygons
        plt.subplot(2, 2, 3)
        
        # Draw the unit circle
        circle = plt.Circle((0, 0), 1, fill=False, color='black', linewidth=2)
        plt.gca().add_patch(circle)
        
        # Draw polygons at specific iterations
        colors = ['red', 'green', 'blue', 'purple', 'orange']
        for idx in [0, 2, 4, 6, min(8, len(results)-1)]:
            n, _ = results[idx]
            points = []
            for i in range(n):
                angle = 2 * math.pi * i / n
                x = math.cos(angle)
                y = math.sin(angle)
                points.append((x, y))
            
            polygon = Polygon(points, fill=False, color=colors[idx % len(colors)], 
                             linewidth=2, label=f'{n} sides')
            plt.gca().add_patch(polygon)
        
        plt.axis('equal')
        plt.xlim(-1.2, 1.2)
        plt.ylim(-1.2, 1.2)
        plt.title('Polygon Approximations')
        plt.legend(loc='upper right')
        plt.grid(True)
        
        # Plot convergence table
        plt.subplot(2, 2, 4)
        plt.axis('off')
        table_text = "Convergence of Polygon Method:\n\n"
        table_text += "Number of Sides | π Approximation | Error\n"
        table_text += "-" * 45 + "\n"
        
        for n, p in results:
            pi_approx = p / 2  # perimeter / 2 for a unit circle
            error = abs(pi_approx - math.pi)
            table_text += f"{n:14d} | {pi_approx:15.12f} | {error:10.12f}\n"
        
        plt.text(0, 0.5, table_text, fontsize=9, fontfamily='monospace')
        
        plt.tight_layout()
        plt.savefig('polygon_pi_approximation.png')
        plt.close()
        print(f"Created visualization 'polygon_pi_approximation.png'")
    
    # Return the final approximation (perimeter / 2 for a unit circle)
    return results[-1][1] / 2


def historical_approximations(visualize=False):
    """
    Calculate π using various historical approximations.
    
    This method implements several historical methods for approximating π
    and compares their accuracy.
    
    Args:
        visualize: Whether to create visualization
        
    Returns:
        Dictionary of historical approximations and their values
    """
    # Dictionary to store the historical approximations
    approximations = {
        "Egyptian (Rhind Papyrus, 1650 BCE)": (16/9)**2,
        "Archimedes (250 BCE)": Decimal('22')/Decimal('7'),
        "Liu Hui (263 CE)": Decimal('3.14159'),
        "Zu Chongzhi (500 CE)": Decimal('355')/Decimal('113'),
        "Aryabhata (499 CE)": Decimal('62832')/Decimal('20000'),
        "Madhava (1400 CE)": Decimal('3.14159265359'),
        "Viète (1593)": Decimal('3.1415926536')
    }
    
    # Create visualization if requested
    if visualize:
        plt.figure(figsize=(12, 8))
        
        # Convert years to numeric for plotting
        year_map = {
            "Egyptian (Rhind Papyrus, 1650 BCE)": -1650,
            "Archimedes (250 BCE)": -250,
            "Liu Hui (263 CE)": 263,
            "Zu Chongzhi (500 CE)": 500,
            "Aryabhata (499 CE)": 499,
            "Madhava (1400 CE)": 1400,
            "Viète (1593)": 1593
        }
        
        years = [year_map[name] for name in approximations.keys()]
        values = [float(value) for value in approximations.values()]
        errors = [abs(float(value) - math.pi) for value in approximations.values()]
        
        # Plot the timeline
        plt.subplot(2, 1, 1)
        plt.plot(years, values, 'o-', markersize=10)
        plt.axhline(y=math.pi, color='r', linestyle='--', label='True π')
        
        # Add value labels
        for i, (year, value) in enumerate(zip(years, values)):
            plt.annotate(f"{value:.8f}", 
                         (year, value), 
                         textcoords="offset points",
                         xytext=(0,10), 
                         ha='center')
        
        plt.xlabel('Year (negative = BCE)')
        plt.ylabel('π Approximation')
        plt.title('Historical Approximations of π')
        plt.grid(True)
        plt.legend()
        
        # Plot the error timeline
        plt.subplot(2, 1, 2)
        plt.plot(years, errors, 'o-', markersize=10, color='green')
        
        # Add error labels
        for i, (year, error) in enumerate(zip(years, errors)):
            plt.annotate(f"{error:.8f}", 
                         (year, error), 
                         textcoords="offset points",
                         xytext=(0,10), 
                         ha='center')
        
        plt.xlabel('Year (negative = BCE)')
        plt.ylabel('Error')
        plt.title('Historical Error in π Approximations')
        plt.yscale('log')
        plt.grid(True)
        
        plt.tight_layout()
        plt.savefig('historical_pi_approximations.png')
        plt.close()
        print(f"Created visualization 'historical_pi_approximations.png'")
    
    return approximations


def buffon_needle_method(needles=1000, visualize=False):
    """
    Estimate π using Buffon's needle problem.
    
    In this method, needles of length L are dropped onto a floor with
    parallel lines spaced distance d apart. If L = d, then π can be
    estimated as 2*n/c where n is the total number of needles and
    c is the number of needles that cross a line.
    
    Args:
        needles: Number of needles to drop
        visualize: Whether to create visualization
        
    Returns:
        Approximation of π
    """
    # Line spacing and needle length (we make them equal)
    d = 1.0
    L = d
    
    crossings = 0
    needle_data = []
    
    for _ in range(needles):
        # Random position of the needle's center
        x = random.uniform(0, 2*d)
        y = random.uniform(0, d)
        
        # Random angle of the needle
        theta = random.uniform(0, math.pi)
        
        # Determine if the needle crosses a line
        # The needle crosses if the center is within L/2*sin(theta) of a line
        distance_to_line = min(y, d - y)
        crosses = distance_to_line < (L/2) * abs(math.sin(theta))
        
        if crosses:
            crossings += 1
        
        # Store the needle data for visualization
        x_end1 = x - (L/2) * math.cos(theta)
        y_end1 = y - (L/2) * math.sin(theta)
        x_end2 = x + (L/2) * math.cos(theta)
        y_end2 = y + (L/2) * math.sin(theta)
        
        needle_data.append((x_end1, y_end1, x_end2, y_end2, crosses))
    
    # Calculate π approximation
    pi_approx = 2 * needles / crossings if crossings > 0 else float('inf')
    
    # Create visualization if requested
    if visualize:
        plt.figure(figsize=(12, 8))
        
        # Draw the parallel lines
        for i in range(-1, 4):
            plt.axhline(y=i*d, color='black', linestyle='-', linewidth=1)
        
        # Plot the needles
        for x1, y1, x2, y2, crosses in needle_data[:min(300, len(needle_data))]:
            color = 'red' if crosses else 'blue'
            plt.plot([x1, x2], [y1, y2], color=color, alpha=0.5)
        
        plt.xlim(0, 2*d)
        plt.ylim(-d/2, 3*d/2)
        plt.title(f"Buffon's Needle Method: π ≈ {pi_approx:.6f} (Error: {abs(pi_approx - math.pi):.6f})")
        plt.xlabel('x')
        plt.ylabel('y')
        
        # Add a legend
        plt.plot([], [], 'r-', label='Crossing Lines')
        plt.plot([], [], 'b-', label='Not Crossing')
        plt.legend()
        
        plt.grid(True)
        plt.savefig('buffon_needle_pi.png')
        plt.close()
        print(f"Created visualization 'buffon_needle_pi.png'")
    
    return pi_approx


def compare_methods(visualize=False):
    """
    Compare different geometric methods for calculating π.
    
    Runs all the implemented methods and compares their accuracy and speed.
    
    Args:
        visualize: Whether to create visualizations
        
    Returns:
        Dictionary of results with method names, values, times, and errors
    """
    print("Comparing geometric methods for π calculation:")
    print("-" * 50)
    
    methods = [
        ("Circle Approximation", lambda: circle_approximation_method(10000, visualize)),
        ("Midpoint Integration", lambda: integration_methods('midpoint', 10000, visualize)),
        ("Trapezoid Integration", lambda: integration_methods('trapezoid', 10000, visualize)),
        ("Simpson Integration", lambda: integration_methods('simpson', 5000, visualize)),
        ("Polygon Perimeter", lambda: polygon_perimeter_method(6, 10, visualize)),
        ("Buffon's Needle", lambda: buffon_needle_method(10000, visualize))
    ]
    
    results = {}
    comparison_data = {'Method': [], 'Value': [], 'Error': [], 'Time (s)': []}
    
    for name, method in methods:
        print(f"\nRunning {name}...")
        start_time = time.time()
        
        try:
            result = method()
            elapsed = time.time() - start_time
            
            # Calculate error
            error = abs(float(result) - math.pi)
            
            print(f"Result: {result}")
            print(f"Error: {error:.12f}")
            print(f"Time: {elapsed:.6f} seconds")
            
            results[name] = {
                'value': result,
                'time': elapsed,
                'error': error
            }
            
            # Add to comparison data
            comparison_data['Method'].append(name)
            comparison_data['Value'].append(float(result))
            comparison_data['Error'].append(error)
            comparison_data['Time (s)'].append(elapsed)
            
        except Exception as e:
            print(f"Error in {name}: {e}")
    
    # Get historical methods too
    try:
        historical = historical_approximations(visualize)
        for name, value in historical.items():
            short_name = name.split('(')[0].strip()
            error = abs(float(value) - math.pi)
            results[name] = {
                'value': value,
                'error': error
            }
    except Exception as e:
        print(f"Error in historical approximations: {e}")
    
    # Create comparison visualization
    if visualize:
        plt.figure(figsize=(14, 10))
        
        # Filter methods for readability
        modern_methods = [m for m in comparison_data['Method'] if 'historical' not in m.lower()]
        indices = range(len(modern_methods))
        
        # Create bar chart for error comparison
        plt.subplot(2, 1, 1)
        errors = [comparison_data['Error'][i] for i in range(len(modern_methods))]
        plt.bar(indices, errors, color='skyblue')
        plt.xticks(indices, modern_methods, rotation=45, ha='right')
        plt.yscale('log')
        plt.ylabel('Error (log scale)')
        plt.title('Error Comparison of π Calculation Methods')
        plt.grid(True, which='both', axis='y')
        
        # Create bar chart for time comparison
        plt.subplot(2, 1, 2)
        times = [comparison_data['Time (s)'][i] for i in range(len(modern_methods))]
        plt.bar(indices, times, color='orange')
        plt.xticks(indices, modern_methods, rotation=45, ha='right')
        plt.ylabel('Execution Time (seconds)')
        plt.title('Performance Comparison of π Calculation Methods')
        plt.grid(True, axis='y')
        
        plt.tight_layout()
        plt.savefig('pi_methods_comparison.png')
        plt.close()
        print(f"Created visualization 'pi_methods_comparison.png'")
    
    return results


def main():
    """
    Main function to run the geometric π calculation demonstrations.
    """
    getcontext().prec = 50  # Set decimal precision
    
    print("Geometric Pi Calculator")
    print("======================")
    
    try:
        # Check if we have command line arguments
        if len(sys.argv) > 1:
            if sys.argv[1] == "--circle":
                points = int(sys.argv[2]) if len(sys.argv) > 2 else 10000
                print(f"Running circle approximation with {points} points...")
                result = circle_approximation_method(points, True)
                print(f"π ≈ {result}")
                
            elif sys.argv[1] == "--integration":
                method = sys.argv[2] if len(sys.argv) > 2 else 'midpoint'
                intervals = int(sys.argv[3]) if len(sys.argv) > 3 else 10000
                print(f"Running {method} integration with {intervals} intervals...")
                result = integration_methods(method, intervals, True)
                print(f"π ≈ {result}")
                
            elif sys.argv[1] == "--polygon":
                sides = int(sys.argv[2]) if len(sys.argv) > 2 else 6
                iterations = int(sys.argv[3]) if len(sys.argv) > 3 else 10
                print(f"Running polygon method with {sides} initial sides and {iterations} iterations...")
                result = polygon_perimeter_method(sides, iterations, True)
                print(f"π ≈ {result}")
                
            elif sys.argv[1] == "--buffon":
                needles = int(sys.argv[2]) if len(sys.argv) > 2 else 10000
                print(f"Running Buffon's needle method with {needles} needles...")
                result = buffon_needle_method(needles, True)
                print(f"π ≈ {result}")
                
            elif sys.argv[1] == "--historical":
                print("Running historical approximations...")
                historical_approximations(True)
                
            elif sys.argv[1] == "--compare":
                visualize = sys.argv[2].lower() == 'true' if len(sys.argv) > 2 else False
                compare_methods(visualize)
                
            else:
                print("Unknown option. Available options:")
                print("  --circle [points]      - Circle point sampling method")
                print("  --integration [method] [intervals] - Numerical integration")
                print("  --polygon [sides] [iterations] - Polygon perimeter method")
                print("  --buffon [needles]     - Buffon's needle method")
                print("  --historical           - Historical approximations")
                print("  --compare [visualize]  - Compare different methods")
        else:
            # Run the comparison by default with visualizations
            print("Running comparison of geometric methods for π calculation...")
            print("(This will generate several visualization images)")
            compare_methods(True)
            
    except KeyboardInterrupt:
        print("\nCalculation interrupted by user.")
    except Exception as e:
        print(f"\nAn error occurred: {e}")


if __name__ == "__main__":
    main()